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ABSTRACT 

Weak lensing measurements are starting to provide statistical maps of the 
distribution of matter in the universe that are increasingly precise and comple¬ 
mentary to cosmic microwave background maps. The probability distribution 
(PDF) provides a powerful tool to test non-Gaussian features in the convergence 
field and to discriminate the different cosmological models. In this letter, we 
present a new PDF space Wiener filter approach to reconstruct the probability 
density function of the convergence from the noisy convergence field. We find 
that for parameters comparable to the CFHT legacy survey, the averaged PDF 
of the convergence in a 3 degree field can be reconstructed with an uncertainty 
of about 10%, even though the pointwise PDF is noise dominated. 

Subject headings: cosmology: observations - cosmology: theory - dark matter - 
gravitational lensing - large-scale structure of universe - methods: N-body simu¬ 
lations 


1. Introduction 

Mapping the mass distribution of matter in the Universe has been a major challenge 
and focus of modern observational cosmology. There are very few direct ways to weigh 
the universe. The only direct procedure to weigh the matter in the universe is by using the 
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deflection of light by gravity. Weak gravitational lensing provides a direct statistical measure 
of the dark matter distribution regardless of the nature and dynamics of both the dark and 
luminous matter intervening between the distant sources and observer. Weak lensing by 
large-scale structure can lead to the shear and magnification of the images of distant faint 
galaxies. While this effect is very small, a large statistical sample can provide a precise 
measurement of averaged quantities. 

Based on the theoretical work done by Gunn (1967), Blandford et al. (1991), Miralda- 
Escude (1991) and Kaiser (1992) performed the Erst calculation of weak lensing by large-scale 
structure, the result of which showed the expected distortion amplitude of weak lensing effect 
is at a level of roughly a few percent in adiabatic cold dark matter models. Kaiser (1992) 
also made early predictions for the power spectrum of the shear and convergence using linear 
perturbation theory. Due to the weakness of the effect, all detections have been statistical in 
nature, primarily in regimes where the signal-to-noise is less than unity. Fortunately, several 
groups have been able to measure this weak gravitational lensing effect (Bacon et al. 2000; 
Refregier et al. 2002; Hoekstra et al. 2002; Van Waerbeke et al. 2002; Jarvis et al. 2003; 
Brown et al. 2003; Hamana et al. 2003) in recent years. 

In the standard model of cosmology, fluctuations start off small, symmetric and Gaus¬ 
sian. Even in some non-Gaussian models like topological defects, initial fluctuations are 
still symmetric: positive and negative fluctuations occur with equal probability (Pen et al. 
1994). As fluctuations grow by gravitational instability, this symmetry can no longer be 
maintained - over densities can be arbitrarily large, while under dense regions can never 
have less than zero mass. This leads to a skewness in the distribution of matter fluctua¬ 
tions. Pen et al. (2003) have measured the first detection of dark matter skewness from the 
VIRMOS-DESCART survey. They find the lensing skewness can be detected to be for a 
compensated Gaussian filter on scales of 5.37. Using the skewness of the convergence for 
simulated weak lensing, Zhang et al. (2003) present the first optimal analysis for the Canada- 
France-Hawaii-Tclescope(hereafter CFHT) Legacy Survey in the presence of noise due to the 
randomly oriented intrinsic cllipticity of source galaxies. They show that a compensated 
Gaussian filter on a scale of 2.5 arc minutes can optimizes cosmological constraint with 
AQ m /Q m ~ 10%, which is significantly better than other window function that have been 
considered in the literature. While skewness has already been measured at very high statis¬ 
tical significance (Pen et al. 2003), the measurement has not resulted in a strong constraint 
on the total matter density Q m . The data has so far been limited by sample variance and 
analysis techniques. Currently ongoing surveys, such as CFHT Legacy Survey, will provide 
more than an order of magnitude improvement in the statistics. 


Besides the skewness of the convergence, the probability density distribution (hereafter 
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PDF) provides a powerful tool to test non-Gaussian features of the convergence field and 
to discriminate the different cosmological models (Jain et al. 2000). In real data of weak 
lensing, the noise due to the intrinsic ellipticity of background galaxies overwhelms the 
PDF of convergence. So far there is no good method to extract the underlying PDF of 
the convergence from the observed noisy version. In this letter, we present a new Wiener 
reconstruction approach to reconstruct the PDF distribution of noise-free convergence held 
from the convergence held added with a Gaussian noise held. 

The outline of the paper is as follows. In §2, we introduce the strategy of map con¬ 
struction of weak lensing from simulations. In §3, we describe the PDF Wiener method of 
reconstructing the distribution of the convergence, while we present the results and summa¬ 
rize our conclusions in §4 and §5 respectively. 


2. Map Construction of Weak Lensing by Simulations 

A full non-Gaussian statistical description of the distribution of matter is very complex. 
It would require the joint n-point distribution function of matter, which is difficult to describe 
or compute. A drastic simplification occurs if we only consider the one point distribution 
of the density held. In gravitational lensing, one only measures the projected dark matter 
density. We then smooth the held with some appropriate hlter, and consider the one point 
PDF of the smoothed density held. Of particular interest is the deviation of this one point 
PDF from Gaussianity. 

A separate question is the degree of information loss implied by this assumption. This 
can be quantified in the context of errors on cosmological parameters in a parameterized 
cosmological model. In a future paper, we will present the accuracy with respect to a set 
of parameters. In this letter, we will simply address the challenge of reconstructing the 
underlying PDF from a noisy measurement. In a controlled simulation, we compare the 
reconstructed PDF to the underlying noise-free PDF. 


2.1. N-Body Simulation 

We ran an N-body simulation with a cosmological model of WMAP to generate con¬ 
vergence maps. The power spectra for given parameters were generated using CMBFAST 
(Seljak & Zaldarriaga 1996) and these tabulated functions were used to generate initial con¬ 
ditions. The power spectra were normalized to be consistent with the earlier two point 
analysis from this data set (Van Waerbeke et al. 2002). We ran the simulations using a par- 



4 


allel, Particle-mesh N-body code PMFAST (Merz et al. 2004) at 1856 3 mesh resolution using 
928 3 particles on an 4 node quad processor Itanium Beowulf cluster at CITA. This publicly 
available N-body code features very fast execution and negligible memory overhead: only 
the positions and velocity of particles are stored, requiring six numbers. One more integer is 
used for every particle to store linked lists. All other arrays are small and recycled. Parallel 
execution efficiency is achieved by decomposing the gravitational force into a long range and 
short range component. The long range force is solved on a coarse grid, which requires very 
little computing or communications resources. The fine grid is only stored in patches, and 
is used to compute the short range force. The highly optimized vendor IPP FFT library is 
used. 

Output times were determined by the appropriate tiling of the light cone volume with 
joined co-moving boxes from z ~ 3 to z = 0. We output periodic surface density maps 
at 1856 2 resolution along the 3 independent directions of the cube at each output interval. 
These maps represent the raw output for the run and are used to generate convergence 
maps in the thin lens and Born approximations by stacking the images with the appropriate 
weights through the comoving volume contained in the past light cone. 

The simulation started at an initial redshift Zi = 80, and ran for 285 steps in equal 
expansion factor ratios with box size L = 100/r -1 Mpc comoving. We adopted a Hubble 
constant h = 0.71, and a flat cosmological model with Q m + Ha = 1 was used. We used the 
WMAP parameters of Q m = 0.27 (Bennett et al. 2003). The power spectrum normalization 
<7 8 was chosen as 0.84. 


2.2. Construction of Simulated Convergence Maps 


The convergence k is the projection of the matter over-density 
6 weighted by the lensing geometry and source galaxy distribution. 


5 along the line of sight 
It can be expressed as 


k(0,Xs)= f W(x)S(x,r(x)0))dx, ( 1 ) 

Jo 

where, x is the comoving distance in unit of c/H 0 , and H 0 = 100 h km/s/Mpc. The weight 
function W(x) is 

W{x) = ^ m g(x)( 1 + z) ( 2 ) 

determined by the source galaxy distribution function n(z) and the lensing geometry. 

^ N f j / , ^r(x'-x) 

9\X) = r{X> / ax n(x )- r7 — 

J x r(x') 


( 3 ) 
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r(x) is the radial coordinate. r(x) = sinh(x) for open, r(x) = X f° r flat and r(x) = sin(%) 
for closed geometry of Universe. In our computations, we adopt a flat geometry. n(z) = 
n (x)dx/dz is normalized such that J 0 °° n{z)dz = 1. For the CFHT Legacy Survey, we adopt 
n(z) = - r ^ + (^) Q ex P(~(~) /3 ) with a = 2 and f3 = 1.2 and the source redshift parameter 
z 0 = 0.44, which peaks at z p = 1.582 0 , respectively. The mean redshift is z = 2.1 z 0 and the 
median redshift is Zh = 1.92o(Van Waerbeke et al. 2002). The source redshift distribution 
n(z) adopted here is the same as that for VIRMOS. 

During each simulation we store 2D projections of 5 through the 3D box at every light 
crossing time through the box along all x, y and z directions. Our 2D surface density sectional 
maps are stored on 1856 2 grids. After the simulation, we stack sectional maps separated by 
a width of the simulation box, randomly choosing the center of each section and randomly 
rotating and flipping each section. The periodic boundary condition guarantees that there 
is no discontinuities between any two adjacent boxes. We then add these sections with the 
weights given by W(z) onto a map of constant angular size, which is generally determined by 
the maximum projection redshift. To minimize the repetition of the same structures in the 
projection, we alternatively choose the sectional maps of x, y, z directions during the stacking. 
Using different random seeds for the alignments and rotations, we make 40 maps for each 
cosmological model. Since the galaxy distribution peaks at z ~ 1, the peak contribution of 
lensing comes from z ~ 0.5 due to the lensing geometry term. Thus the maximum projection 
redshift 2 ~ 2 is sufficient for the lensing analysis. So we project the Q m = 0.27 simulations 
to 2 = 2 and obtain 40 maps each with angular width 0 K = 3 degrees. At the very highest 
redshift bins, the field of view is larger than the simulation volume. In those slices we used 
periodic boundary conditions to fill in the field. Their lensing contributions are negligible. 
Fig.(l) shows a noise-free k map in the N-body simulation of a WMAP cosmology with a 
width of 3 degrees and 2048 2 pixels, and the scale is in units of k. 

We then simulate the CFHT Legacy Survey by adding noise to these clean maps. The 
noise k maps have a pixel-pixel variance a 2 N = (e 2 )/2/(W pixel ) = 0.74 2 , which is the noise 
variance in each pixel before wiener Liter. Here (e 2 ) = 0.47 2 is the total noise estimated 
in the VIRMOS-DESCART survey and here we take it as what would be expected by the 
CFHT Legacy Survey. It includes the dispersion of the galaxy intrinsic cllipticity, PSF 
correction noise and photon shot noise. (At pixel ) is the mean number of galaxies in each pixel. 
For VIRMOS, the number density of observed galaxies n g ~ 26/arcmin 2 , then (V pixel ) = 
n g[(6 K /l')/N] 2 , where N = 2048 is the number of grids by which we store 2D maps and the 
Leld of view 0 K is in units of arc min. The factor of 2 arises from the fact that the shear held 
has two degrees of freedom ( 71 , 72 ), where the dehnition of (e 2 ) sums over both. We use this 
as our best guess for the CFHT Legacy Survey noise. We smooth both of the noise-free and 
noisy convergence Lelds by a Wiener Liter Fffi, so %(i) = f k(x')Fw(x — x'). The Fourier 



6 



Fig. 1.— A noise-free k map in the N-body simulation of a WMAP cosmology with a width 
of 3 degrees and 2048 2 pixels, and the scale is in units of k. 
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transform of the Wiener filter is given by 

F ^ k) = wmhk) (4) 

where P s (k) is the power spectrum of the k field and P N (k ) is the power spectrum of the 
noise. The broad Wiener filter smoothes the image, and reduces the variance of the noise 
per pixel. This can be revealed by the dotted line in Fig. (4), where the variance of pure 
Gaussian noise field smoothed with the Wiener filter is = 0.0084. The Wiener filter has 

an effective area of {cfnI^n/) 2 pixels. Of course, the pixels are highly correlated now. 

The real k, field is quite non-Gaussian, and a Wiener filter might not be the best. An 
alternative decomposition based on wavelets is discussed by Pen (1999). 

Fig. (2) and (3) show the the noise-free and noisy k maps by adding simulated Gaussian 
noise field smoothed by Wiener filter respectively. It is apparent that even the minimum 
noise Wiener maps are very noisy, and where noise dominates salient features in the maps. 
The k field is equivalent to the E-mode of the shear field ( 7 ), and a Wiener filter allows a 
direct mapping of the shear to filtered n field (Padmanabhan et al. 2003). It is well known 
that shear and kappa have equivalent noise properties (Seljak 1998). This can be proven as 
follows: white and independent noise in qx and 72 , when Fourier transformed, is still white. 
The decomposition into E and B modes, i.e. the component of the shear in Fourier space 
that is parallel and perpendicular to the wave vector, is just a rotation of two uncorrelated 
noise variables. This gives two new uncorrelated noise variables of equal variance. The E 
mode and B mode will both have white noise again, which is statistically identical to the 
noise in 71 and 72 (differing by only a rotation). The E-mode is identical to kappa, except 
for one constant mode. Therefore kappa and shear can have their noise treated identically. 
For real surveys, the naive Wiener filter costs as the number of galaxies to the third power, 
and could be computationally very expensive, but fast O(N) algorithms have already been 
demonstrated (Pen 2003). 


3. Wiener Reconstructing of the Convergence PDF 

We apply Wiener filtering in two completely different contexts. In the previous sections, 
the Wiener filter was used to generate a minimum noise variance map. The next step is the 
tabulation of the PDF from a 9 square degree simulated survey field. This PDF is necessarily 
a noisy realization of the ensemble expectation value. 

In addition to the sample variance noise, the observations themselves have noise. This 
is noise is added pointwise to the intrinsic lensing surface density. The PDF of the sum of 




Fig. 2.— The same as Fig.l, but for the noise-free k map smoothed by a Wiener filter. 
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Fig. 3.— The same as Fig.2, but for noisy k map by adding a simulated Gaussian noise field 
smoothed with the same Wiener filter. 
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signal plus noise is the convolution of the two PDF’s. We model the noise PDF as a Gaussian 
with zero mean and specified standard deviation. We further assume that the noise standard 
deviation is known. 

In the ensemble average, one could in principle recover the lensing PDF from the noisy 
PDF by a deconvolution. In practice, of course, this is very difficult because the small scale 
noise is amplified by the deconvolution to the extent that it completely overwhelms the 
signal. Lee (2002) showed a heuristic scheme for noisy Sunyaev-Zeldovich measurements, 
which regularizes this deconvolution by applying a Wiener filter on the PDF of peaks. Lee 
extracted the “power spectrum” of the signal and noise in the PDF using Fonrier transforms. 
Fourier transforms are good descriptions for stationary processes, but clearly a PDF is not 
a stationary process as a function of k. 

Stationarity is not a requirement for Wiener filtering. One can specify the full covari¬ 
ance matrix. However, one must be careful to formulate the question. What is meant by 
covariance? What is one approximating by a Gaussian process? 

We first consider the literal interpretation. We define the ensemble average PDF of the 
40 maps as the underlying PDF. The actual PDF of each map is a noisy realization. The 
difference of each map PDF from the ensemble average is the noise free signal variance. 

Next, we consider a single map PDF with noise added. If the process of adding noise were 
an exact convolution, this would be a reversible process for a non-zero kernel, in particular 
for a Gaussian. To assess the actual contribution of noise, we take the noise-free PDF P s , 
and convolve by the ideal noise kernel. Then we take the noisy map, and measure its actual 
noisy PDF P s+N . The difference of the noisy map’s PDF from the ideal convolved PDF is a 
measure of the noise’s contribution on the irreversibility of the convolution. We thus define 
the noise covariance matrix as the covariance of the difference between the ideal convolved 
PDF and the noisy PDF of each map. The Wiener filtered deconvolution can then be readily 
applied to the observed noisy PDF, as we will demonstrate below. 

One might wonder how circular this procedure really is. Our ultimate goal is to find 
a minimum noise measurement of the ensemble PDF from a noisy measurement of a single 
held PDF. In order to implement this, we had assumed prior knowledge of an ensemble of 
fields, as well as knowledge of the ensemble average k PDF and its covariance! At first sight 
this seems to imply much more prior knowledge than we set out to solve for. In practice, the 
problem is not as serious as it might sound. It needs to be shown that the actual assumed 
PDF for the signal is not important, or that it can be derived from the observed data alone. 


We thus test not only the literal interpretation of a Wiener reconstructed PDF, but also 
of a more liberal version, we have assumed the signal PDF to be Gaussian. For a Gaussian, 
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only the variance needs to be known or measured from the data, which can already be done to 
good accuracy from existing lensing data. The quality of the reconstruction will necessarily 
be worse, but we will show from actual simulations that good results can nevertheless be 
achieved. 


The PDFs of pure k held, pure noise and k held added by noise are P s (k), P n (k) 
and P s+n (k) respectively. The number of maps is IV = 40, i.e. n — 1,2, • • ■ N. Normally, 
PDFs are normalized to be unit, i.e. P(k)(Ik = 1. From here on we only consider 

binned PDFs, choosing M equally spaced bins separated by An. Our notation becomes 
Ki = K- + iAn, Pi = P(hii)An, and each PDF can also carry an appropriate superscript. 
The number of bins for all PDFs is M = 19, i.e. i,j — 1, 2, • • • M. 

In the form of matrix, C s is the signal covariance and C N is a noise covariance. The 
signal covariance is expressed as 


Cf,=< \pf -Pf}[Pf --Pf]>, 


( 5 ) 


where <> is the average over 40 maps, and P s (n) is the average value of the k held PDF for 
40 maps. The noise PDF, which will also serve as the noise convolution kernel g N (n — n'), 
is a Gaussian distribution 



1 


^2707 jv 


exp 


(Ki ~ Kj ) 2 

2 o% 


( 6 ) 


where On is the standard deviation of the noise after the wiener filter. We can also think of 
this kernel as a noise matrix N 


Nij = g N (ni - Kj) Ak. (7) 

For a noisy k held, the noise deviation of the PDF relative to an invertible convolution can 
be written as for each map 

ApN = pS+N _ NijP S, (8) 

where P s+N (ni) is the PDF of k held with added noise. Thus we have the noise covariance 
matrix 

C*j =< [AP t N - AP™} [APf - A Pf ] > . (9) 


Using the convolution theorem, we have 


r + OO 


P s+N = / P S (K,')g N (K- K')dK, 


M 


p, s+N = Y, N ‘r p i’ 

3 = 1 


( 10 ) 


( 11 ) 
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which reduces to in the form of matrix 


P S+N = NP 5 , 

(12) 

so we have 


pS _ N-lpS+iV 

(13) 

The Wiener filter is defined as, 


w = C 5 (C 5 + C Ar )- 1 , 

(14) 


where C N = N^C^N” 1 . 


All covariances are specified as deviations from a known underlying model. Given a the 
PDF of a noisy map, we then must first subtract the ensemble average noisy PDF, and apply 
the Wiener reconstruction to difference. After a applying a Wiener filtered deconvolution, 
we add back the noise free ensemble average PDF. The figure of merit is how well the 
reconstruction finds the noise-free PDF of the original map (not the ensemble average). 

The reconstruction is thus decomposed into two parts: an underlying mean noisy PDF 
has a known pre-convolved (i.e. deconvolved) PDF. The difference between the map realiza¬ 
tion and this ensemble average is then deconvolved with a Wiener filter. 

Expressed as equations, the deviation of reconstructed P 5 17 is 

AP 5 - = C S N[NC S N + C^AP 5 ^ (15) 

where 

ap s+n = P s+N _ p . P = NP . (16) 

Thus we can obtain the reconstructed P^ 

P^ = AP^ + P 5 , (17) 


and its error 


°ls 


=< (P 


P 5 ) 2 >, 


w 


where <> is the mean over 40 maps, and P s is the PDF of the k field. 


(18) 


A key assumption was that we actually knew the ensemble average PDF. The difference 

between the observed PDF and the ensemble average is a smaller amplitude function, and 

can be deconvolved using an appropriate Wiener filter. In practice, one might not know the 

underlying PDF, and subtract the wrong function. To test this possibility, we also used a 

~s 

Gaussian mean distribution for P in equation (16). In the reconstructed PDF we again 
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used the exact deconvolution of the Gaussian (which is of course also a Gaussian) P . In 
this case, we are applying the exact deconvolution to a lesser component of the observed 
PDF, and test how important it is to know the true underlying ensemble average PDF. One 
could in principle also use the signal covariance field of a Gaussian field in Equation (5). It 
does not appear as crucial a qualitative difference, and was not tested. 


4. The reconstructed probability density function (PDF) and the cumulative 

distribution function (CDF) 


The PDF of a convergence map is defined as 


P( K ) 


ti(Ak) 

N 2 Ak 


(19) 


where N is the number of grids by which 2D maps are stored and ti(k) is the number of 
the within unit interval Ak of k. Thus its corresponding cumulative distribution function 
(CDF) can be expressed as 

C(> k) = f P{k )(1k. (20) 

J AC 

where P(k) is normalize to be unit in our analysis. Our purpose is to extract the PDF 
of k from noisy k map and compare it with that of pure convergence field, thus we can 
test the method mentioned in last section and apply it to the analysis of real week lensing 
data. Using the matrix method, we reconstruct the PDF of convergence from the pure, 
noisy k map and pure Gaussian noise smoothed by a Wiener filter. Fig.(4) plot the PDFs 
of the noise-free k field (solid line), noisy k held (dashed line) and pure Gaussian noise held 
(dotted line) respectively. The corresponding error bars are obtained from averaging over 
the 40 respective maps. The empty square dots represent the reconstructed PDF of the k 
held using the simulated pure k held as the initial signal held. We see that the reconstructed 
PDF of k is in good agreement with that of simulated pure k held, while there exists slight 
deviation for smaller negative and larger positive k. Fig. (5) shows their corresponding CDFs. 


In the reconstruction process, we used the PDF of the simulated pure k held as a 
initial signal held. However, in the real observation of weak lensing, the PDF of the pure 
convergence is unobservable. Thus we employ an approximation that the PDF P s (kj) of 
simulated convergence is replaced by a Gaussian distribution with mean standard deviation 
(J K over 40 simulated k helds, because we can infer the mean standard deviation of k from the 
observational data of weak lensing. Fig. (6) and (7) plot the reconstructed PDFs and their 
corresponding CDFs using a Gaussian held with a mean standard deviation of 40 simulated 
pure k maps as an initial signal held, while Figs. (8) shows the relative errors of the PDF 
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for the reconstructed hi field using an initial signal field of a simulated pure hi field (solid 
line) and a Gaussian field (dashed line) respectively. Apart from in smallest negative and 
largest positive convergence bins, the relative errors for the two initial PDF distributions are 
similar and about below 10%. Therefore, using a Gaussian held PDF with a mean standard 
deviation of 40 simulated pure n maps as an initial signal held is a good approximation in 
the reconstruction of convergence PDF using this Wiener reconstruction approach. 


5. Conclusions and Discussions 

In this letter, we ran a high-resolution N-body simulation of a WMAP cosmology to 
study the reconstruction of the convergence PDF from noisy weak lensing data. We added 
noise due to intrinsic ellipticity of background faint source galaxies to the simulated hi helds 
and smoothed it using a Wiener hlter. From the noisy simulated hi held, we make the 
reconstruction of pdf of convergence by means of the PDF Wiener method described above. 
We find the reconstructed PDF of the convergence are in good agreement with that of a 
noise free n map smoothed by a Wiener hlter, and its relative error is below about 10%. 
We can safely apply this reconstructed method to the analysis of real observational data for 
week lensing such as future CFHT Legacy Survey. 
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